Fluctuations of internal energy flow in a vibrated granular gas. 
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The non-equilibrium fluctuations of power flux in a fluidized granular media have been recently 
measured in an experiment [Phys. Rev. Lett. 92, 164301, 2004], which was announced to be a 
verification of the Fluctuation Relation (FR) by Gallavotti and Cohen. An effective temperature 
was also identified and proposed to be a useful probe for such non equilibrium systems. We explain 
these results in terms of a two temperature Poisson process. Within this model, supported by 
independent Molecular Dynamics simulations, power flux fluctuations do not satisfy the FR and 
the nature of the effective temperature is clarified. In the pursue of a hypothetical global quantity 
fulfilling the FR, this points to the need of considering other candidates than the power flux. 

PACS numbers: 02.50.Ey, 05.20.Dd, 81.05.Rm 



Granular gases, i.e. gases of macroscopic grains losing 
part of their kinetic energy during collisions, have been 
the subject of ardent investigation in the last decade Q. 
They display, in experiments as well as in simulations, a 
rich and intriguing phenomenology: non-Gaussian statis- 
tics, breakdown of energy equipartition, spontaneous 
symmetry breaking (clustering, shear, convection, sur- 
face waves, shocks etc) are the most striking features 
of this surprising state of matter (see e.g. [2( and ref- 
erences therein) . In addition, a driven granular gas is an 
ideal testing ground for non-equilibrium statistical me- 
chanics of dissipative stationary states. An inelastic gas 
can be kept in a state of constant average total kinetic 
energy thanks to an external driving mechanism, e.g. by 
rapidly vibrating its container. In such a system, en- 
ergy is flowing from the external source/thermostat into 
the gas and then from the gas into an irreversibly drain- 
ing sink, represented by inelastic collisions. As a conse- 
quence, the fluctuations of global or microscopic physical 
observables, such as the total energy or the velocity of 
particles or all internal currents, do not behave as ex- 
pected in equilibrium statistical mechanics. It is there- 
fore tempting to give an interpretation of the observed 
fluctuations in terms of the few known theoretical results 
in non-equilibrium statistical physics. In a recent exper- 
iment on vibrated granular gases Q it has been argued 
that the statistics of the power flux fulfilled the Fluctua- 
tion Relation (FR) by Gallavotti and Cohen 0,0, even if 
-as emphasized in |3|- a) the FR holds under conditions 
that are not met in the situation under study (in partic- 
ular, microscopic reversibility is required) and b) there is 
no proven relation between the measured power flux and 
the entropy production entering the original FR. 

We show here that these experimental results can be 
explained in terms of a simple Poisson process. We argue 



that the measured "power flux" Q is well reproduced by 
a sum of random and independent energy amounts char- 
acterized by two different typical temperatures. Such a 
quantity does not verify the FR. We nevertheless demon- 
strate that a straight line with slope /3 e // can be ob- 
served when plotting \og[f Q (Q) / f Q {-Q)] vs. Q, (f Q 
being the probability density distribution of Q), in the 
range which is accessible in experiments and simulations. 
In |3] the predictability of Peff remained an open prob- 
lem and 1 / P e f f was interpreted as an effective temper- 
ature. Within the two temperature model, the value of 
(3 e f f is obtained and it further appears that interpreting 
this quantity as an inverse effective temperature is prob- 
lematic. We validate our analysis by a direct confronta- 
tion against Molecular Dynamics (MD) simulations, ob- 
taining very good agreement. 
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FIG. 1: Left: Snapshot of the system considered for MD sim- 
ulations, with the inner region marked by the solid rectan- 
gle. Right: Corresponding vertical profiles of density ($(?/), 
dashed line) and temperature (T(y), solid line). The dotted 
lines mark the bottom and top boundaries of the inner region. 
Here N = 270 and a — 0.9. The mean free path is ~ 5.7d. 

The event driven MD simulations have been performed 
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for a system of N inelastic hard disks with restitution 
coefficient a, diameter d and mass m = 1. The verti- 
cal 2D box of width L x = 48d and height L y = 32d is 
shaken by a sinusoidal vibration with frequency / (pe- 
riod Tb ox = 1//) and amplitude 2.6d. In a collision two 
particles lose a fraction 1 — a 2 of their relative kinetic 
energy while the total momentum is conserved. Colli- 
sions with the elastic walls inject energy and allow the 
system to reach a stationary state. We have checked 
that possible inelastic collisions with the walls hardly af- 
fect the results. Gravity - set to g = — 1.7df 2 in or- 
der to be consistent with the experiment- has a neg- 
ligible influence on the measured quantities. We have 
varied the restitution coefficient from 0.8 up to 0.99 
(glass beads yield on average a w 0.9) and the total 
area coverage from 0.138 (i.e. N = 270) up to 0.32 
(N = 620). In Fig. QJ-left a snapshot of the system 
is shown. During the simulations all the main phys- 
ical observables are statistically stationary. The local 
area coverage field 3>(x,y) and the granular temperature 
field T(x, y) (defined in 2D as the average kinetic en- 
ergy per particle) are almost uniform in the horizontal 
direction, apart from small layers near the side walls. In 
Fig. fright the profiles $(?/) = (1/L X ) j dx<fr(x,y) and 
T(y) = (1/L X ) J dxT(x,y) are shown to be symmetric 
with respect to the bottom and the top of the box. Fol- 
lowing the experimental procedure, we have focused our 
attention on a "window" in the center of the box, fixed in 
the laboratory frame, of width 2L x /5 and height L y /3, 
marked in Fig.GJ-left. Apart from the negligible change of 
potential energy due to gravity, the total kinetic energy 
of the particles inside the window, changes during a time 
r because of two contributions: AK T = Q T — W T where 
Q T is the kinetic energy transported by particles through 
the boundary of the window (summed when going- in and 
subtracted when going-out) and W T is the kinetic en- 
ergy dissipated in inelastic collisions during time r. For 
several values of r we have measured, as in the experi- 
ments, Q T which is related to the kinetic contribution to 
the heat flux (we checked that inclusion of the collisional 
contribution, even if non small [6j, does not change the 
picture). With N = 270 and a = 0.9 the characteristic 
times are the mean free time t co ; w 0.47t;, ox , the diffusion 
time across the window Tdiff = 0.82t{, ox and the mean 
time between two subsequent crossings of particles (from 
outside to inside) T cross w 0.039t{, ox . 

We define the injected power as q T = Q t /t and two rel- 
evant probability density functions (pdfs): /q{Qt) and 
f q (q T )- Fig. EK) shows f q {q T ) for different values of r. 
A direct comparison with Fig. 3 of Ref. 0] suggests a 
fair agreement between simulations of inelastic hard disks 
and the experiment. The pdfs are strongly non-Gaussian 
and asymmetric, becoming narrower as t is increased. 
At small r a strong peak in q T = is visible. More 
interestingly, f q (q T ) at small values of r has two differ- 
ent exponential tails, i.e. f q (q T ) ~ exp(=p/3-trg T ) when 
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FIG. 2: a) pdfs of injected power f q (q T ) from MD simulations 
for different values of r = (1, 2, 4, 8, 16, 32) x r m i„ with r m i„ = 
0.015rbox- Here iV = 270 and a — 0.9. The distributions are 
shifted vertically for clarity. The dashed lines put in evidence 
the exponential tails of the pdf at r = T m tn- b) plot of 
(1/t) log[/ 9 (q T )//g(— Qt)] vs. q T from MD simulations (same 
parameters as above) at large values of r. The solid curve 
is a linear fit (with slope Peff) of the data at r = 128r m in. 
The dashed line has a slope f3 gra n = i/T gran . In the inset the 
same graph is shown for small values of r = (1, 2, 4, 8) x r m i„ 
(from bottom to top). 



q T — > ±oo with f3- > /?+. The peak and the expo- 
nential tails at small r are observed also in the exper- 
iment (see Fig. 3 of and in similar simulations 0. In 
Fig. Hb) we display \og[f q {q T )/f q (-q T )]/T vs. q T , which 
is equivalent to the graph of 7r(g r ) — 7r(— q T ) vs. q T where 
ir(q T ) = log[/ Q (rg T )]/T. When r -> oo, n(q T ) -> 11(g), 
i.e. the large deviation function associated to /q(Q t )- 
Under a number of hypothesis, the Gallavotti-Cohen 
Fluctuation Theorem |j| states that for the entropy pro- 
duction fj (defined in a dynamical system as the phase 
space contraction rate) n(er) — II(— a) = a. From Fig. [21 
it appears that at large values of t, Tr(q T ) — 7t(— q T ) is 
linear with a r-independent slope /J e // ^ 1. We have 
measured /3 e / / with various choices of the restitution co- 
efficient a and of the covered area fraction finding sim- 
ilar results. Ref. 0] reports f3 e ffT gran ~ 0.25 where 
Tgran is the mean granular temperature in the observa- 
tion window. Similar values are measured in our MD 
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simulations. At area fraction 13.8% and a = 0.9 we have 
(3 e ffT gran ?s 0.23. At fixed a and increasing area frac- 
tion, PeffTg ran slightly increases, as in the experiment. 
As a — > 1 the slope e ff decreases. At a = 1 (without 
gravity and external driving) the distribution of Q T is 
symmetrical and /3 e // = 0, indicating that l//3 e // is not 
a physically relevant temperature concept. Interestingly, 
it appears that e ff is a non hydrodynamic quantity: 
different systems may show the same density and tem- 
perature profiles, with very different values of /3 e // 01 . 

We now adopt a coarse-grained description of the ex- 
periment which is able to entirely capture the observed 
phenomenology. The measured flow of energy is given by 
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where n_ (n+) is the number of particles leaving (en- 
tering) the window during the time r, and vf_ (vf + ) are 
the squared moduli of their velocities. In order to ana- 
lyze the statistics of Q T we take n_ and n + as Poisson- 
distributed random variables with average cot, where ui 
corresponds to the inverse of the crossing time T cross . 
In doing so we neglect correlations among particles en- 
tering or leaving successively the central region. A key 
point, supported by direct observation in the numerical 
experiment, lies in the assumption that the velocities v,-+ 
and come from populations with different tempera- 
tures T_|_ and T_ respectively. Indeed, compared with 
the population entering the central region, those parti- 
cles that leave it have suffered on average more inelastic 
collisions, so that T_ < T + . Finally we assume Gaussian 
velocity pdfs || . Within such a framework, the distribu- 
tion Jq{Q t ) of Q T can be studied analytically. Here it is 
enough to recall that ^ X)"=i v 1 ■> m D dimensions, if each 
component of Vj is independently Gaussian-distributed 
with zero mean and variance T, is a stochastic vari- 
able x with a distribution Xn,T(x) = fi/T,Dn/2{x), where 
fa,u{x) is the Gamma distribution, and whose generating 

function reads Xn,r{z) = (1 — Tz)~ Dn ^ 2 0. It is then 
straightforward to obtain the generating function of Q T 
in the form /q(z) = cxp[Tfi(z)] with 



fj,(z) = Lu[-2 + (l-T + z) 



-D/2 



(1 + T.z) 



-D/2 
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We observe that Jq{z) has two poles in z = ±1/T± and 
two branch cuts on the real axis for z > 1/7+ and z < 
— l/TL. From we immediately obtain the cumulants 
of f Q (Qr) through the formula (Q n ) c = r^/^O). 

For t — > oo the large deviation theory states that 
Iq{Qt) ~ exp(rn(Q T /r)) and H(q) can be obtained 
from n{z) through a Lcgendre transform, i.e. H(q) = 
max(/i(z) — qz). The study of the singularities of n{z) 
reveals the behavior of the large deviation function IL(q) 



for q — > ±oo. It can be seen that 

n(g)~-^- («-»«>), n(g) 



(g--oo). (3) 



We emphasize however that it is almost impossible to 
appreciate these tails in simulations and in experiments, 
since the statistics for large values of qand r is very poor. 

A Gallavotti-Cohen-type relation 0, , e.g. 11(g) — 
II(— q) = (3q for any q and an arbitrary value of (3 would 
imply /i(z) = /i(/3 — z). One can see that such a f3 does 
not exist, i.e. the fluctuations of Q T do not satisfy a 
Gallavotti-Cohen-like relation. The observed linearity of 
the graph log[f q (q T )/ f q (-q T )]/T = ir(q T ) - n(-<k) vs - 
q T can be explained by the following observation at 
large values of r it is extremely difficult, in simulations as 
well as in experiments, to reach large values of q, while for 
small q, ir(q)— 7r(— q) « 2ir' (0)q+o(q 3 ), i.e. a straight line 
with a slope /3 e // = 2ir'(0) is likely to be observed. It has 
been already shown |l0| that in dissipative systems devia- 
tions from the FR can be hidden by insufficient statistics 
at high values of q. The knowledge of n(z) is useful to 
predict this slope. At large t, tt'(O) w IT(0) = -z*(0) 
where z*(q) is the value of z for which n{z) — qz is ex- 
tremal. This gives 



7-11 T + 
e ff = 2 tt=- with 7 = — — ; S = — 

1 1 7 + -f 5 T- T_ 2 
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When 7 = 1 (i.e. if a = 1) e ff =0. As a decreases, 
7 increases, since the collisions dissipate more energy, 
and e ffT_ grows, reaches a maximum and subsequently 



decreases asymptotically toward as 
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phasize that e ff does not depend upon ui. We have 
compared these predictions with the numerical and ex- 
perimental results, measuring the temperatures T + and 
T_ in the simulation. To bypass the problem of tempera- 
ture anisotropy (discussed below) , we have used the hor- 
izontal component of the temperature, obtaining (with 
a = 0.9 and area fraction 13.8%) T+ « Mid 2 / 2 and 
T_ « 91d 2 / 2 , i.e. 7 = 1.55 and, from Eq Q in D = 2, 
e ff = 0.00193 in very good agreement with the mea- 
sured value 0.0022. It should be noted that the tempera- 
ture ratio can also be approximated by (a 2 )~ m where m 
is the average number of collisions undergone by a parti- 
cle between the moments of entering and leaving the ob- 
servation region. In our numerical simulations (as well as 
in the experiment) m « 2, which corresponds to 7 « 1.5. 

What happens for small values of r? We note that 
/q(z) has the form cxp(r/i(z)) for any value of t 
and not only for large r. Therefore at small r one 
can expand the exponential, obtaining /q(z) ~ 1 + 
cut (-2 + (1 - T + z)- D l 2 + (1 + T_z)- D / 2 ). This imme- 
diately leads to an analytical expression for /q(Q t ) = 
const x S(Q T ) + X\,t+{Qt) + Xi,T-(—Qt), which fairly 
accounts for the strong peak which is observed in the 
experiment and in the simulations, and predicts expo- 
nential tails for /q(Q t ): Xi.t{x) oc x d ^ 2 ~ x exp(— x/T) 
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so that (3+ = 1 /T+ and j3- = 1 /T_ . This suggests an ex- 
perimental test of this theoretical approach: the measure 
at small values of t of the slopes of the exponential tails of 
Iq(Qt) should coincide with a direct measure of T + and 
T_. However, we point out that the values of f3 + and /3_ 
obtained by fitting the tails in the hard disks simulation, 
using values as small as t = 0.00015t£, o:e yield estimates of 
T + and Tl which are larger (by a factor ~ 1.6) than those 
found by a direct measure. This disagreement brings the 
limits of such a simple two-temperature picture to the 
fore. In the simulation and in the original experiment the 
measured injected energy is indeed the sum of several dif- 
ferent contributions, namely Q T « Q% x +Q%* +Qf> +QW 
where is the kinetic energy transported by the i com- 
ponent of the velocity by particles crossing the bound- 
ary through a wall perpendicular to direction j. Two 
main differences with the simplified interpretation given 
above arise: a) there are two couples of temperatures, 
i.e. T%,T* as well asT|,T! jUj; b) the diagonal con- 
tributions Qp are sums of squares of velocities whose 
distribution is not a Gaussian but is ~ v exp(— v 2 /T), 
since the probability of crossing is biased by the velocity 
itself. The calculation of fp(Q T ) is still feasible, with 
qualitatively similar results 12]. 




-2e+06 2e+06 4e+06 

[d¥] 

FIG. 3: Color online: Pdf of transverse (xy) injected 
power f q (q T ) for low values of r = (1,2,4,8) x T m ;„ with 
Tmin = 0.00015 Tbox (resp. circle, squares, diamond and tri- 
angles) from MD simulations with N = 270 and a = 0.9. 
The distributions are shifted vertically for clarity. The (red) 
solid lines show the solution of the "two temperatures" model. 
The dashed lines indicate exponential tails ~ exp(^(3±rq T ). 
Left inset: same plot for a large value r = 6400r m i„. Right 
inset: plot of (1/t) log[f q (Q T )/fq(— Qt)] vs. q = Q/t from 
MD simulations (same parameters as above) at large values of 
r = (1, 2, 4, 8) x T min with r' min = 3200r m i„, together with a 
dashed line of slope /3 e f f predicted by Eq. Q and the dotted 
line of slope gra n = 1/T gran . 

Here we focus on the proposed two-temperature model, 
showing its ability to explain the statistics of internal 
energy currents. To this extent we have repeated the 
MD simulations discussed above, but measuring only 
the transversal kinetic energy current Q* y through the 
bottom boundary of the central region. For new mea- 



surements, shown in Fig. we have also improved 
the time resolution of the measure, using a minimum 
t = 0.00015T cross . The coefficients of the exponen- 
tial tails of f q (q T ) at low t are in perfect agreement 
with the predicted t/T± and the slope of the graph 
log[ f q (q T ) / f q (—q T )]/T vs. q T , shown in the inset, is ac- 
curately recovered by Eq. |J2J with D = 1. More remark- 
ably, the two-temperature model (solved by numerical 
inversion of /) fully accounts for the pdf f q {q T ) at any 
value of r, as evidenced in the same figure. 

In conclusion we have implemented MD simulations 
of inelastic hard disks fluidized in a vibrated box. We 
have studied the statistics of the power flux through a 
closed perimeter -which effectively separates a high tem- 
perature region from a low temperature one- thereby ac- 
curately reproducing the experimental measures of 0]. 
Inspired by these results, we have put forward a simple 
two-temperature Poissonian model, which entirely repro- 
duces the phenomenology of the system but for which the 
fluctuations of the energy flow do not satisfy a Gallavotti- 
Cohen relation. Moreover, our approach puts forward 
a kinetic interpretation for the measured slope /3 e //, a 
question previously left open. We conclude here that 
there are serious evidence against the power flux as a 
potential candidate for extending the FR to a granular 
gas, a system for which the validity of such a relation is 
therefore still open. 
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